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We present numerical calculations of the parameters Ci/, Ce and Ck associated with 
^ the common closures for turbulence production, dissipation and diffusion. In the case 

►^ . of homogeneous and isotropic turbulence, these parameters are expected to be statisti- 

(^ ' cally scale-invariant within the inertial subrange. In order to scrutinise this conjecture, 

OO . we utilised a generalisation of the Germano filtering formalism, which is applicable to 

^^ ' compressible flows as well. The filtering of data obtained from three-dimensional direct 

numerical simulations of forced isotropic turbulence with Mach numbers in the range 
'^ 0.1 ... 1 then yielded values of the closure parameters associated with different length 
f~^ , scales. The results indicate that the closure parameters are nearly universal for subsonic 

or moderately transonic flows, although the resolution of 432'^ grid cells in our simu- 
lations is not quite sufficient to clearly establish scale invariance. In addition, it was 
found that the customary assumption of a kinetic Prandtl number of about unity for the 
gradient-diffusion closure is flawed due to the misalignment between turbulent flux and 
the gradient of the turbulence energy. Nevertheless, sound correlation can be achieved 
^ , if the flux magnitude rather than the flux vector is locally matched. This conclusion is 

particularly useful for the family of subgrid scale models based on the turbulence energy 
equation. Furthermore, the parameter of production Ci, was computed in the fashion 
^\f ' of dynamical procedures. Thereby, superior agreement between modelled and explicitly 

evaluated turbulence stresses in comparison to the eddy-viscosity closure with constant 
Ci, was verified. 



1. Introduction 

The notion of self-similarity plays an important role in the numerical computation of 
turbulent flows, particularly, in large-eddy simulations (LES; seCiSagaut, 2001). The key 
idea is that the large-scale features of the flow are shaped by the action of external forces 
and the imposed boundary conditions, whereas structures on smaller scales are mostly 
determined by the non-linear transfer of kinetic energy through a cascade of vortices. 
The dynamics on the largest scales, which constitute the energy- containing range, must 
be explicitly computed in a LES. On the other hand, turbulent velocity fluctuations on 
smaller scales, which are contained in the so-called inertial subrange, are considered to be 
statistically self-similar. Usually, this property of inertial-range turbulence is expressed 
in terms of scale-invariant probability distribution functions associated with structural 
properties of the flow. It becomes beautifully manifest in the k~^/^ power law of Kol- 
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mogorov's theory. In the context of closures for second or higher order moments of the 
velocity field (and, possibly, other quantities), the associated parameters are expected 
to become asymptotically scale-invariant, once the cutoff length scale sufficiently pene- 
trates the inertial subrange. If this supposition was validated, it would bear important 
consequences on the treatment of the unresolv ed fraction of turbulence in LES, i. e., the 
problem of subgrid scale (SGS) modelling fsee iMeneveau fc: KatJ.l20Qfl) . 

Inspired by the turbulence energy equation model fsee ISagautl l200l[ §4.3.2), we at- 
tempted to explicitly evaluate the rate of production and dissipation, respectively, as well 
as turbulent diffusion in consecutive bands of wave numbers, using data obtained from 
numerical simulations of forced isotropic turbulence. The dissipation of kinetic energy on 
the smallest resolved scales in these simulations is entirely due to the implicit viscosity 
produced by the numerical scheme. The physical viscosity of the fluid is negligible. From 
structural properties of the flow on larger scales, which appear to be fairly independent of 
the details of dissipation, we were able to compute the parameters Ci, and C^, stemming 
from the eddy-viscosity and the gradient-diffusion closure, respectively. For the rate of 
dissipation, we invoked the hypothesis of local equilibrium and computed the parameter 
Ce associated with the simple dimensional closure. Furthermore, the correlation between 
modelled and explicitly calculated quantities was investigated. The numerical simula- 
tions which provided the fundamental data are briefly outlined in 321 The underlying 
formalism of filtered quantities is introduced in ^ The following 21 33 ^-nd ^ cover 
the discussion of various results for the closure parameters and comparisons between the 
corresponding modelled and actual quantities. In the conclusion, the question whether 
we can indeed speak of self-similarity properties inferred form simulation data is posed. 
The analysis presented in this paper points towards a positive answer. However, more 
computing power in order to achieve higher resolution would be desirable. This becomes 
particularly clear from samples of the turbulence energy spectrum function discussed in 
121 In particular, the spectra display the notorious bottleneck effect, which is a deviation 
from inertial-range scaling on length scales clos e to the cutoff spoils the full visibility of 
self-similarity to some extent (see iDobler. Haugen. Yousef fc Brandenburg. 2003) . 



2. Forced Isotropic Turbulence 

In order to produce numerical realisations of isotropic turbulence, we computed the 
flow of a compressible fluid set into motion by a random force field in a cubic domain 
subject to periodic boundary conditions. In the case of statistically isotropic and station- 
ary forcing, the outcome is an almost perfect realisation of the mathematical paradigm 
of homogeneous and isotropic turbulence. Only the imposed periodicity in Cartesian 
coordinates induces an inherent anisotropy on large scales. The dynamics of the fluid is 
generally determined by the following set of hydrodynamical conservation laws: 

(2.1) 

P + pf^ + T—(T^k, (2.2) 

OXk 

(2.3) 

Energy is injected into the fluid through a mechanical force field f{x, t) called the driv- 
ing force and ultimately dissipated by the microscopic viscosity v. The corresponding 
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dissipation tensor atk is proportional to the local rate of strain of the velocity field: 

a,k = 2pvS*t, = 2pv (s,k - \d5,t\ , (2.4) 

where Sik — 2{'Vi,k ^'Vk,i) and d = Vi,i. Since mostly velocity fluctuations on the smallest 
length scales contribute to the strain, microscopic viscous dissipation becomes negligible 
on larger scales and the fluid dynamics is dominated by non-linear turbulent interactions. 
At sufficiently high resolution, the cutoff arising from the discretisation in numerical sim- 
ulations falls into this very range of scales. In consequence, either a subgrid-scale model 
has to be employed in order to account for the energy transfer from numerically resolved 
towards unresolved length scales, or a dissipative finite-volume scheme is applied, which 
properly smoothes the flow on the smallest resolved scales. We cho se the latter approach 
and adopted the piece-wise parabolic method (PPM) proposed bv IColella fc Woodward! 
( (l984) for the solution of the hydrodynamical equations. The idea of employing the 
PPM for direct numerical simulations of turbulent fl ows was particularly advocated by 
ISvtine. Porter. Woodward. Hodson fc Winklen (|200Cf), following systematic convergence 
tests and comparisons to conventional methods with explicit treatment of the viscosity 
term. From the numerical point of view, the treatment of the hydrodynamical equations 
with the PPM corresponds to the limit of zero physical viscosity. In essence, the com- 
pressible Euler equations are solved, with the velocity fluctuations on scales smaller than 
the grid resolution A being damped out by numerical dissipation. 

The driving force f{x, t) is composed in spectral space, using a three-dimensional 
genera lisation of the scalar Ornstein- Uhlenbeck process, as proposed by Es waran fc Pop g 
l|1988|) . The evolution of the Fourier transform f{k,t) is given by the following Langevin- 
type stochastic differential equation: 

1 /? 
dfik,t)^~f{k,t)^ + FoJ2(^^^) 5{k~k,ir^)Pcik)-dWu (2.5) 

jlni 

The second term on the right hand side accounts for a random diffusion process, which 
is constructed from a three-component Wiener process W*. The distribution of each 
component is normal with zero mean and variance dt. The wave vectors kjim are dual 
to the position vectors of the cells in the numerical discretisation of the fundamental do- 
main. The symmetric tensor P({k) is defined by the linear combination of the projection 
operators perpendicular and parallel to the wave vector. The components of P^(fc) can 
be expressed as 

(P,,)c(fc) - CPi^ik) + (1 - C)PII (fc) = C% + (1 - 20^, (2.6) 

where the spectral weight C determines whether the resulting force field in physical space 
is purely solenoidal, dilatational or a combination of both. The variance cr^(fe) specifies 
the spectrum of the force field. We use a quadratic function, which confines the modes 
of the force to a narrow interval of wavenumbers, k G [0, 2fco]. The wave number feg 
determines the integral length scale of the flow, L = l-KJkQ. 

The root mean square of the speciflc driving force is determined by the characteristic 
magnitude Fq and the weight C: 

/r.„s = Y.(hi^^) ■ himit)) ^ (1 - 2C + 3e)Fl (2.7) 

jlrn 

Since Fq has the physical dimension of acceleration, it can be expressed as the char- 
acteristic velocity V of the flow divided by the integral time scale, which is given by 
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NAt 


1.0 


0.084 


2.5 


5815 


1.0 


0.42 


3.0 8.0 


6343 


0.75 


0.66 


5.0 10.0 


6351 


0.2 


1.39 


5.0 10.0 


3356 



Table 1. List of simulation parameters: spectral weight ( of the driving force, characteristic 
Mach number V/cq, onset of decay t^/T, end of simulation tf/T, number of time steps NAt- 



the auto-correlation time T of the driving force 1)2. 5|l . Setting T = L/V, we have 
-fo = V/T = LV^, and, starting with a homogeneous fluid at rest, the flow is devel- 
oping towards a fully turbulent steady state within about two integral time scales. 

As our work was motivated by the problem of turbulent burning processes in ther- 
monuclear supernovae, we applied the equation of state for a degenerate electron gas in 
combination with non-degenerate nucl ei. This form of matter occurs in compact stel - 
lar remnants called white dwarfs (see iHillebrandt fc Niemeveil l200Ci iReineckel l200l|) . 
However, once quantities are scaled in terms of characteristic parameters, no major dif- 
ferences to turbulence in, say, an ideal gas are found. The parameters chosen for several 
particular simulations are listed in tabled In the simulations with (^ = 1, turbulence 
was produced by purely solenoidal forcing or "stirring" . We chose two different values 
of the characteristic Mach number V/cq , where cq is the speed of sound in the initially 
homogeneous fluid at rest. For the lower value, the flow is completely subsonic, whereas 
for the higher value, it is marginally transonic. Furthermore, we run two simulations 
with partially dilatational forcing, i. e., the force fleld has compressive components. In 
both cases, the velocity of the flow locally exceeds the speed of sound and shocklets are 
formed. For all simulations, 432^ grid cells were used and a = X/L = 3, where X is the 
linear size of the fundamental domain. The evolution was computed over an elapse of 
several integral time scales, such that data dumps corresponding to a statistically sta- 
tionary state with energy injection balanced by dissipation were obtained. Finally, decay 
phases were initiated at a certain time td by inhibiting the random increments in the 
evolution of the stochastic force field. Also included in table ^ are the total durations 
if of the simulations in units of the corresponding integral time scales and the required 
total number of time steps. 

Samples of the normalised energy spectrum function E{ak,t) = {2TT/aLV'^)E{k,t) 
at representative stages in the production regime, in statistical equilibrium and in the 
advanced decay regime are plotted in figure ^ i^~'^) for the simulation with C, = 0.75. 
Due to the discrete nature of the numerical data, the kinetic energy of modes within 
certain wave number bins was summed up. Also shown are the longitudinal and a 
transversal fractions of the energy spectrum functions corresponding to, respectively, the 
compressible and incompressible components of the velocity field. Figure Q] (d-f) shows 
the corresponding plots of the compensated solenoidal spectrum function 
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where (e) = (T/V'^){e) is the normalised mean rate of dissipation. In the inertial sub- 
range, ^s{t) ~ C is expected. C is known as the Kolmogorov constant and considered 
to be more or less universal. From the spectra shown in figure^ it becomes clear that 
a resolution of A^ = 432 is just at the brink where an inertial subrange begins to take 
shape. Even in the quasi-equilibrium state, there is merely a narrow window of wave 
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Figure 1. Plots of the transversal and longitudinal parts of the turbulence energy spectrum 
function E{k, t) (a-c: top panels) as well as the compensated transversal spectrum function 
^^(fe, t) (d-f: bottom panels) at different times in a simulation of forced isotropic turbulence. 
The characteristic Mach number is V/co = 0.66 and the spectral weight of the driving force 
C = 0.75. 

numbers in the vicinity of fc = (L/27r)fc = 3.0, in which nearly Kolmogorov scaling with 
C « 1.7 is found. In fact, this value of the Kolmogorov cons tant is c omfortably within 
the range of numerical results reported in the literature (see lYeung fc Zhou . 19) . Both 
in figure n (e) and^ (f), there is a pronounced maximum of the compensated spectrum 
function at fc « 15 corresponding to a flattening of the energy spectrum in comparison to 
the Kolmogorov law (figures ^b an dHlc). This so-call ed bottleneck effect was observed 
in many numerical simulations (see lDobler et a^J . l2003|) . At time i = t/T = 1.5, on the 
other hand, small-scale features have only partially developed and the spectrum does not 
show a bump at higher wave numbers (figures ^^^ and^d). We will not further discuss 
the bottleneck effect here, but take it as a genuine feature of our simulations. 



3. Hierarchical filtering 

In this section we shall elaborate the formal concepts of filters and physical fields 
smoothed over the associated length scales. The formalism was first introduced by 
|Germano| l|1992j) in the framework of incompressible flows. In the following, mass- 
weighted filtered quantities are defined, which are applicable to compressible flows as 
well. Remarkably, one of the key results obtained by Germano, the so-called Germano 
identity, carries over to the compressible case without modification. In order to avoid 
confusion, we will subsequently write the quantities given by the hypothetical exact solu- 
tion of the dynamical equations in the continuum limit, i. e., for infinite resolution, with 
the superscript oo set on top of the corresponding symbol, whereas the standard symbols 
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refer to numerically computable quantities. Notwithstanding the deplorable lack of any 
proof of existence, the former shall be called the ideal quantitiesf. For instance, v{x,t) 
is identified with the ideal velocity field, i. e., the exact solution of the Navier-Stokes or 
Euler equation, while it is understood that v{x,t) is an approximation to the velocity 
field, say, an interpolating function to some finite-volume solution. 

Since we are mostly concerned with statistically stationary homogeneous turbulence, it 
is convenient to introduce an infinite series of homogeneous and time-independent filters 
( )„. The kernel of the n-th filter is denoted as Gn{x)- The corresponding /i/ter operation 
on a dynamical quantity q{x, t) is defined by 



g^ 



(")(a;,f)= d?x'G,,{x')%{x\t) (3.1) 



and symbolically written as (?'■"•' = {q)n- In the limit n ^ oo, the identity operator is 
obtained. Here we shall assume that ( )„ is a Gaussian filter with a characteristic length 
A„ and wave number fc„ = 7r/A„. The corresponding kernel is given by fsee lPooelbOOCl 
§13.2) 



6 \^/' / 6\x[ 



In the following, a series of Gaussian filters with A„ — L/n, n e N, is used, where L is 
the integral length scale of the flow. 

The velocity field t;["l(a;, t) at the n-th filter level is defined by a Favre mass-weighted 
filter operation: 

^'"' (^' ^) = („J ,, / d'^' Gr,{x')^{x', t)v{x', t) (3.3) 

or, for brevity, v'"' = (pv)n/(p)n- Filtering twice, we set 

^M[«]^M(„) ^ (pM^H)„ = {(pv)^)„, (3.4) 

where p(™)(") — {{p)rn)n- If two filters of characteristic length scales A,„ and A„ » A,„ 
are operating in succession, we have 

^H[n]^(m)(„) ^ (pM^H),^ ^ (^^)^^ ^ t;["V(»). (3.5) 

The validity of this approximation becomes immed iately clear from the product of the 
kernels of the Gaussian filters in spectral space fsee lPoDell200Cll §13.2): 



Gni{k)Gn{k) = exp 



fc'(A^ + A2 



24 



exp 



24 



G„(fc). (3.6) 



Because the convolution with the filter kernel in physical space corresponds to the mul- 
tiplication of the Fourier transforms in spectral space, it follows that {{q)m)n — {q)n if 
A„ > A„. 

The definition of the filtered velocity at the level n H3.3|l is particularly useful, because 
the equation of motion for the fiuid on scales I > A„ is given by a quasi- Navier-Stokes 

f According to Plato's Cave Analogy, the objects we experience in the world are related to 
their ideals like the silhouettes of outside things to someone trapped in a dark cave. In the 
same way, one can think of a numerically simulated flow as being merely a vague image of the 
corresponding physical flow or exact mathematical solution of the equations of motion. 



equation: 

The specific force /^ stirring the fluid is assumed to vary on the largest scales only. This 
is why the force field is only marginally affected by the filtering operation, and one can 
set fl o::i fi for n ^ 1. In comparison to the physical Navier-Stokes equation, 

00 

U occtt O 000000 Ul r ^ /o o\ 

-^P^i + ^ — pviVk^-j. ^pjt + -^ — (Tik, (3.8) 

at oxk oxi oxk 

the viscous stress tensor Oik — 2 p vS*f^ for a fluid of microscopic viscosity v is enhanced 
by the turbulence stress tensor associated with the n-th fllterf . The latter is defined by 

rt'=-Cpv^.)n + P^-->vYV:\ (3.9) 

In particular, if A„ is large compared to the length scale of viscous dissipation, then rj^ 

dominates over (J^\ and ( I3.7|) becomes a quasi-Euler equation, in which the microscopic 
viscosity v does not appear at all. 

For two consecutive filter levels, say, n and n — 1, there is a simple algebraic relation 
between the corresponding stress tensors: 



"^ifc 



;^"""'' = (r!i:')„-i + r["--ii(«!"Ul"^), (3.10) 



where 



,0000, , , ,oc 00, 



\n]\n-l] ,,000000 {{pVi)n)n-l{{pVk)n)n-l 

Tik =- {{pViVk)n)n-l-\ ^ = 

{(p)n)n-l (3.11) 

- ((pS?5°fe)„)„_i +p(")("-l)«["l["-%i"l["-ll. 



and 



r[-"-il (-!"', 4"') = -(p(").!"14"')-i + -^;r^(p(")t'!"l)„-i(p(")4"')-i (3.12) 

is the intermediate stress tensor at the filter level n associated with the vel ocity field 
filtere d at the level n — 1 . The above relation was originally formulated by iGermand 
l|1992fl for incompressible flows. The relation also applies to arbitrary filter levels, m and 
n, say. In the Hmit A„ :^ A,„, the contribution from {tI^ )„ becomes negligible and 

r|,"l c, 4"""1 c, r["'"l(z;j"l, 4"'), (3.13) 

i. e., the turbulence stress associated with the scale A„ is not sensitive to the flow 
structure on much smaller scales. The asymptotic limit of the Germano relation is 
especially useful for the numerical evaluation of turbule nce stresses fr om numerical data. 
Utilising the consistent Germano decomposition (see ISaeautl I2OOII §3.3.2), the turbu- 
lence energy associated with length scales I < A„ is given by the difference between the 
filtered and the resolved kinetic ene rgy at the filter level n. This particular notion of 
turbulence energy was introduced bv lGermand l|1992|) and is called the generalised turbu- 
lence energy. At any filter level, the energy is readily defined by the trace of turbulence 

f Here the opposite sign as customary in most of the literature is used in order to make Tik 
a proper stress tensor, which enters the right-hand side of the equation for the resolved energy 
with positive sign. 
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stress tensor: 

Ki^]=p{n)k[n]=_lrlf. (3.14) 

This definition avoids several difficulties stemming from the mathematical properties of 
common filter operations, if the turbulence energy is defined in terms of velocity fluctua- 
tions relative to the filtered velocity field in the alternative Leonard's decomposition (see 
ISagautLl200ll §3.3.1). Contracting the compressible Germano identity H3.1U|I . we obtain 
the relation 

where 

is the intermediate turbulence energy, i. e., the kinetic energy of modes which are con- 
centrated in the spectral band [7r/A„_i,7r/A„]. As opposed to the spectral filter, there 
are no sharp boundaries between adjacent bands associated with Gaussian filters. Nev- 
ertheless, the notion of turbulence energy as proposed by Germano is well-defined and 
unambiguously associates some energy contents with each band of wave numbers. 

The scale separation of the energy budget by means of filtering on a characteristic 
scale A„ also yields a conservation law for the turbulence energy K^"^. The result is 
completely analogous to the dynamic al equation for SGS turbulence energy and it also 
entails various closure problems (see iGermand Il992t ISaeautL 120011 S3. 3 . 2). A dopting 
the standard closures for production, diffusion and dissipation fsee lSagauul200ll §4.3.2), 
we obtain the following approximate equation for the time evolution of K^"^': 

°* ''' ^ - ., . (3.17) 

where ^y^ = ^ -I- t>["l • V is the Lagrangian time derivative with respect to the filtered 
velocity field, 5*^. = ^{dkv'P' + div'j^') — |d'"li5ifc are the components of the trace-free 

rate-of-strain tensor and rfl"! = div^ is the divergence of the velocity field filtered at the 
level n. Now the question of self-similarity boils down to the scaling-behaviour of the 
parameters Cl^ , Ce and Ck associated with the different filters in the hierarchy. 



4. Kinetic energy transfer 

The rate of transfer of kinetic energy from velocity fluctuations on length scales I > A„ 
towards those on smaller scales is given by the contraction of the turbulence stress tensor 
and the rate of strain tensor at the level of the n-th filter: 

nH^rl;:!^!:!. (4.1) 

In] In] 

Note that S^j^ is the symmetrised derivate of the filtered velocity field, i. e., S'lj, = 

2{dkv" + 9ivj." ). Since the filter operation involves mass- weighing, differentiation and 
filtering do not commute. The most common closure for the rate of energy transfer is 
the eddy -viscosity closure for the trace-free part of the turbulence stress tensor: 

4l"l=2p(").(")5;J"l, (4.2) 



where r*^,'"^ = r|^' — ^T^Sik = rj^ + ^K^^^^Sik, and z^^"^ is the turbulent viscosity of the 
fluid at the length scale A„. Since viscosity can be expressed as the product of a length 
scale and a characteristic velocity, a customary hypothesis identifies j^*^"' = Ci" A„\/fcH. 
Hence, 

jj(n) _ 2^NrfN ^ p(")ci")A„\/fcN|S'*["lp. (4.3) 

o 

If turbulence were self-similar within a certain range of length scales, then one would ex- 

(n) 

pect Cy to be scale invariant, when it was averaged over appropriate regions of the flow. 
In order to test this proposition a priori, we filtered simulation data in the intermediate 
range of scales between the grid resolution and the integral scale. For the explicit evalu- 
ation of the turbulence stress tensor, an enhanced viscosity approximation was applied. 
This means that the ideal velocity field in definition (|3.9|) is replaced by the numerically 
computed field, which is smooth on scales / < A due to numerical discretisation. Equiva- 
lently, one can think of the smoothness being caused by an associated numerical viscosity, 
which enhances the physical viscosity of the fluid. By the same line of reasoning as in the 
case of two filters with A„ ^ A„i, it follows that (g)„ — {{q)cff)n = {q)n if A„ ::§> Aog. 
Moreover, 3N : Ajv+i < Acff < Ajy. According to the asymptotic equation (|3.13|) . rj^' 
can be approximated by r^j, ' . This, in turn, implies 

4"lc.-(p.,.,)„ + p(").["l.["l (4.4) 

for a filter of significantly larger characteristic length than the the numerical scale Acff. 
Here Vi is the PPM solution of the quasi-Euler equation H2.2|l . 

However, the actual range of inertial scales in the case of the simulation discussed in 
121 is rather marginal. From the turbulence energy spectrum shown in figure ^ one can 
see that approximate Kolmogorov scaling is found for dimensionless wave numbers in the 
narrow range 2 ^ fc ^ 5. The dimensionless wave number associated with a Gaussian 
filter of characteristic length scale A„ is fc„ — (7r/A„)(L/27r) = L/2A„. Thus, only filters 
with 0.1 ^ An/L ^ 0.25 are more or less suitable for calculating Ci" from H4.3|) . if the 
the turbulence stress tensor is substituted by the enhanced viscosity approximation (|4.4|) . 

(n) 

Notwithstanding these tight constraints, mean values of Ci, were calculated from a 
sample of data sets using several different filters. The results are plotted as functions of 
the characteristic filter wave number in figure|21 Although, as a consequence of the rather 
limited resolution of the simulations, pronounced self-similarity is not apparent, basic 
trends can be discerned. In the case of purely subsonic turbulence with the characteristic 
Mach number V/cq « 0.084, the plotted values seem to indicate a maximum of (Ci" ) 
near the wave number k = 2.0, which marks the upper bound of the energy-containing 
subrange. Towards higher wave numbers, i. e., for filters of smaller characteristic length, 
{Cv) decreases and eventually flattens in the vicinity of maximum dissipation at the 
wave number fcp « 13.5. For the other simulations, a similar behaviour emerges, but 
there is seemingly a trend toward smaller production parameters for increasing Mach 
number and partially dilatational forcing. In any case, we suggest to adopt {C^,) ~ 0.06 
as a good value for fully developed turbulence. This choice is further supported by the 
selection of numerical values listed in table |21 which were obtained with the particular 
filter of characteristic wave number fcg = 3. Near this wave number, the compensated 
energy spectra exhibit local minima. 

A different notion of self-similarity originates from the correlation of turbulence stresses 
modelled upon the eddy-viscosity closure in consecutive bands of wave numbers. Actu- 
ally, so-called dynamical procedures for the computation of subgrid-scale closure param- 
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Figure 2. Numerically evaluated closure parameters (Cu), Ct and Ck as functions of the 
normalised wave number fc„ — L/2A„ for Gaussian filters of characteristic length A„ . 
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Table 2. Mean production, dissipation and diffusion parameters calculated with a Gaussian 
filter of characteristic length Ag — L/6 = 24A. The values in the first three columns specify 
the characteristic Mach number V/cq, the spectral weight ^ of the driving force and the instant 
of time t/T for the chosen data sets extracted from different numerical simulations. 



eters, in particular, C^, make use of this kind of similarity. Originally, dynamical formu- 
lations of the well-known Smagorinsky model and the SGS turbulence energy equation 
mo del utilised the proposition that turbulence stresses at different filter levels are siniilar 
CseelOermano. Pio mclU. M oin fc Ca bof. 1991; E rlebacher. Hussaini. Sneziale fc ZangLll992t 
IPiomelli fc LiuL Il995i: .Meneveau fc Kata. |2(}00). The ratio of the filter length scales in 
this approach is usually set equal to a factor of two. Thus, using the notation introduced 
in 13 the similarity hypothesis amounts to C'"-* = (C'-^"-')„. Contrary to the statistical 
scale invariance of C'"' investigated earlier in this section, this is a local relation, which 



presumes significant correlation between the turbulence stresses Tj). and (r^j, )„ 



(n) 



^(2n)\ 
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On the other hand, iLiu. Meneveau fc Kata 1)1994 ) found evidence for a more pro- 
nounced correlation between r|j,"' and the intermediate stress tensor rj^,"'" . The latter 
is also known as the Leonard's stress tensor in the context of SGS models. They pro- 
cessed data obtained from velocity measurements in round jet experiments with filters 
of varying characteristic length scale in order to evaluate the turbulence stresses and the 
rate of energy transfer across certain wave numbers. Making use of the data produced i n 
our numerical simulations, we re- investigated the relation reported bv lLiu et al\ l|l994J) . 
To this end, the function Ci "■'" (a;,t) was determined from the intermediated stress 
tensor rlfT associated with the range of scales A2„ ^ / £ A„ = 2A2„: 

The quantities in the above expression can be evaluated by filtering the data at two levels 
n and 2n, in between the energy containing and the dissipation range. The similarity 
closure for the turbulence stress tJ^" associated with the length scale A2„ is then given 

by 

^*J2„]^[2„] ^ C(2n,n)^(2n)^^^ypi;;] 1^* [2n] |2^ (4 g) 

In order to validate the modelled rate of production according to (|4.6|) . once more the 
turbulence stress tensor on the left-hand side was evaluated by means of the enhanced 
viscosity approximation. 

The deviations of both the statistical and the similarity closure from the explicitly 
evaluated rate of turbulence energy production is illustrated in figure |31 First of all, 
figure ini (a) shows a contour plot of r* '^'5'|J, which was computed with a Gaussian filter 
of characteristic length Ae — L/Q from a 2D section of the flow in the simulation with 
purely solenoidal forcing and characteristic Mach number V/cq « 0.42 at time t ~ 2.5. 
The corresponding turbulent-viscosity closure with the production parameter set equal to 
(Ci ) « 0.0615 is plotted in figure|31(c). Although the overall agreement is satisfactory, 
regions of pronounced production are evidently not reproduced. On the other hand, the 
similarity parameter Ci ' was computed from (|4.5|l . The corresponding contour plot is 
shown in figure O (b) . Obviously, Ci is negative in several regions of the flow. The 
resulting negative turbulent visc osity is usually interpreted as hackscattering of energy 
from smaller to larger scales fsee ISaeautl l200ll §4.4). In fact, energy transfer upwards 
through the cascade is well known from turbulence theory and can be accounted for by 
localised closures. The outcome of inserting C\j ' into the turbulent-viscosity closure 
for r* [^1 S\f} is shown in flgure (d) . The deviations from explicitly evaluated rate of 
production are clearly less than in the case of the closure with the constant statistical 
parameter. Although some residual errors persist, the localised closure performs much 
better in regions of large turbulence stresses. Moreover, it gives a fair approximation to 
backscattering. 



5. The rate of dissipation 

The rate of dissipation at the level of the n-ih fllter, 

oo 

e(")=z.(|5*|2)„, (5.1) 

specifies the rate of conversion of kinetic energy into internal energy due to microscopic 
viscous dissipation. The filter operation corresponds to local averaging of the rate of 
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FlGURE 3. The rate of turbulence production evaluated from filtered data of the simulation 
with purely solenoidal forcing and V/cq « 0.42 in the plane z; = at time t = 2.5 (a: left top 
panel) and the corresponding deviations of the statistical and the localised turbulent-viscosity 
closure, respectively (c, d: bottom panels). Also shown is the localised production parameter 
inferred from a similarity hypothesis (b: right top panel). 



dissipation over length scales I < A„. Regardless of the filter length A„, the rate of 
dissipation e*^") is largely determined by velocity fluctuations on the smallest dynamical 
scales near the microscopic Kolmogorov length ijk- For this reason, there is no way of 
explicitly evaluating e'-"-' from under-resolved numerical data. A simple estimate can be 
made if the resolved flow is more or less in stochastic equilibrium. Taking the global 
average of the turbulence energy equation H3.17|l . the diffusion term cancels out, and the 
mean time derivative is negligible in the case of steady turbulence. Hence, 



_i_n(«)^c. (.(")+ AN). 
The contribution A*^"^ stems from pressure dilatation: 



(5.2) 



(5.3) 



Equation 1)5. 2|l expresses the well known equilibrium hypothesis for the balance between 
turbulence production and dissipation in the quasi-stationary regime for the dynamics 
on length scales I > A„. The validity of this hypothesis has recently been confirmed for 
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incompressible turbulenc e by directly computing the spectral energy transf er fmiction 
from numerical data fsee lKaneda. Ishihara. Yokokawa. Itakura fc UnoLl2flfl.'Il . 

The most common closure for the rate of dissipation is the dimensional expression 

6(")=C,(") \^ / / . (5.4) 

In general, the parameter Ce varies in space and time. However, in the case of sta- 
tistically stationary homogeneous turbulence, the closure hypothesis can be simplified 

(n) 

further by assuming a constant value of Ce ■ In order to fulfil statistical equilibrium, 
we therefore set 

Ci"-' = ^r 3/2\ ■ (5-5) 

Results calculated from simulation data are shown in figureEl It appears that the graphs 

(n) 

of Ce flatten near the transition to the range of wave numbers dominated by dissipation. 
In particular, for the simulations with lower Mach number, Ce ~ 0.48. Even for the 
higher Mach numbers, the corresponding values of the dissipation para meter are not 
much different. This falls in place very well with the results discussed bv lKaneda et all 
l|2003|) . A further sample of values for the dissipation parameter is listed in tabic |2] In 

(n) 

most cases, Ce is slightly less than 0.5 for the filter levels 6 < n < 12, which highlights 
the robustness of this parameter. 



6. Turbulent diffusion 

The turbulent flux of kinetic energy is given by 

4"' = ^^+Mr\ (6.1) 

where the contracted third-order moments (TOM) t^^j^ and the pressure-diffusion flux 
fj,^"^ are defined by 



T, 



(n) ^/Jp JP !^\ /°S????!?!^\ o^["L.["l 

iik 



T{Vi,Vi,Vk) = -{pViViVk)n-^r^k Vi + {P^i^i/nVk, (6.2) 

and 

Aife = T(P,^fe) = _(P^.)„ + F(")4"1. (6.3) 

In the well known gradient-diffusion closure, the flux is set equal to the mean product of 
kinetic diffusivity and the gradient of the turbulence energy: 

fI"! =p(")C^")A„v/A^Vfc["l. (6.4) 

(n) 

This closure immediately reveals a consistency problem. The parameter C^ is over- 
determined, as there are three flux components. For a well deflned solution, the gradient 
of /c^"! being aligned with the flux vector F^"'' is an indispensable precondition. If this 
were the case, a least-square approach could be employed in order to calculate the closure 
parameter: 

, , (j^M . VA:["1) 

C(") = ^ ^ ' . (6.5) 

(p(«)A„\/fcR|VfcN|2) 

In the particular case of the simulation with partially dilatational forcing {(, — 0.66), 



14 



W. Schmidt, W. Hillebrandt and J. C. Niemeyer 



^Q'l.F^^F?Y^/fiy 



^^^^T^y-^T^TV'^^.^.^'^. 




10ViL,C^ [k[°'kMk^)VVp,v' 



0. 1 £ 5 

'A 




Figure 4. The turbulent flux magnitude computed with a Gaussian fiher of characteristic 
length Ag — L/6 from data of the simulation with (" = 0.75 and V/co ~ 0.66 at time i — 5.0. 
Shown are the contours of the actual flux magnitude \F^^'\ (a: left top panel), the projection 
of the modelled flux as defined by equation 16.51 unto the actual flux (b: right top panel), and 
contours of differently modelled fluxes (c, d: bottom panels) . 



the outcome of this ansatz was matched against the explicitly evaluated flux F^ ' for a 
2D section of the filtered flow at time i = 4.0. The filtering length is Ag = L/6. Equa- 
tion Ht).5|l yields Ck ~ 0.0358 from the full 3D data set, which is of the same order of 
magnitude as the corresponding mean parameter of production, C^ w 0.0574. From 
this, a turbulent kinetic Prandtl number close to one seems to be supported. However, 
the overall agreement between the modelled and the explicitly evaluated turbulent flux is 
actually very poor, as is revealed by a comparison of the figures0|(a) and^(b). The flux 
magnitude is underestimated by almost a decade, according to the different scales of ei- 
ther plots. This discrepancy definitely indicates that the a priori assumption of turbulent 
transport in the direction of the turbulence energy gradient is not fulfilled. In particular, 
those points where the turbulence energy gradient is oriented nearly perpendicular to 
the turbulent flux vector appear as dark wiggling ribbons in figure 0] (b) . 

In consequence, we forsook the approach outlined so far and determined the diffusivity 
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parameter by matching the flux magnitudes rather than the flux vectors: 

Ci"> . ^S3 ^ (6^6) 



Implementing this equation, C^ ~ 0.390 was obtained from the aforementioned sim- 
ulation data. In fact, this value is larger by about an order of a magnitude than the 
corresponding value computed with (|6.5|) . Contours of the corresponding modelled flux, 

p'^^'^Ck Agvfcl^l Vfc[^l|, are plotted in figure 0] (c). The remarkably good correlation 
to I-F' 'l is evident. Even surfaces at which the flux vanishes, as a result of pressure- 
diffusion cancelling the TOM contributions, are well reproduced. Therefore, we conclude 
that the gradient-diffusion closure is a fair statistical description of turbulent transport, 
which correctly accounts for the magnitude but not for the local direction of transport. 
Numerically calculated values of Ck for different simulations are listed in table 13 It 
appears that there is a trend towards stronger diffusion for higher Mach numbers. Fig- 
ure [21 shows the variation of Ck with the smoothing length scale A„. As one can see, 
the parameter of diffusion is almost scale-invariant for n € 2, . . . , 10, which corresponds 
to the nearly inertial subrange. 

In addition, we tested whether turbulent diffusivity is correlated with the turbulent 
viscosity, since k^") — cr^j^i'^") with (Tkin ~ 1 is a commonly used hypothesis. To this end, 
the parameter Ci" was computed locally, according to the prescription in 21 Then the 
diffusivity was set equal to the viscosity for positive values of Ci" . On the other hand, 
Cy was set identically zero in regions of negative production, because concentration 
of turbulence energy due to negative diffusivity has to be inhibited. The numerical 
evaluation demonstrated that the resulting diffusive flux is typically too small by an 
order of a magnitude. An exemplary 2D contour section is shown in figure 01 (d). The 
assumption of a kinetic Prandtl number close to unity is thus clearly invalidated. 



7. Conclusion 

We have investigated parameters associated with the common closures for the turbu- 
lence energy budget equation, using data from direct numerical simulations of compress- 
ible isotropic turbulence driven by stochastic forcing. The statistical self-similarity of 
isotropic turbulence, which becomes manifest in the scale-invariance of averaged closure 
parameters within the inertial subrange, just begins to emerge at the resolution of 432 
grid cells per dimension. Nevertheless, it was feasible to infer the estimates (C^) w 0.06, 
Ce ~ 0.5 and C^ ~ 0.4 for the parameters of turbulence production, dissipation and 
diffusion, respectively. These estimates appear to be robust for Mach numbers < 1 and 
predominantly solenoidal forcing. Significant deviations from these values are likely to 
emerge from supersonic turbulence, although we are not able to make definite statements 
on the basis of our data. Whereas the results for the parameters of production and dis- 
sipation, respectively, agree well with various values cited in the literature, we propose a 
substantially enhanced parameter of diffusion. The would-be smallness of C^, as given 
by the standard hypothesis C„ ~ C^, can be attributed to the mismatch arising from 
the usually presumed vector alignment in the gradient-diffusion closure. We emphasise 
that there is no a priori reason for any tight relationship between the turbulence viscos- 
ity (specifying the quasi-local transport of energy from larger to smaller scales) and the 
turbulence diffusivity (specifying the non-local redistribution of energy on the smaller 
scales). 
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Furthermore, we compared the modelled turbulence stresses for a localised eddy- 
viscosity closure to the explicitly evaluated stresses and found improved correlation, in 
contrast to the the closure with a constant statistical parameter. This result supports the 
validity of dynamical procedures for the computation of Ci, in SGS models , espec ially, 
when applying the subgrid scale model proposed bv iKim. Menon fc Mongial I 1999|) . 

For a more stringent analysis of the self-similar regime in numerical simulations, sig- 
nificantly higher resolution is called for. Indeed, extremely high-resolution data have 
been computed recently (sec Kaneda et at . 2003') . These data could be exploited for 
the evaluation of structural properties over a wide range of length scales by means of 
filtering. 
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